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Low  Observability  Path  Planning  for  an  Unmanned  Air 
Vehicle  Using  Mixed  Integer  Linear  Programming 

Atif  Chaudhry,  Kathy  Misovec,  and  Raffaello  D’ Andrea 


Detection  of  an  Unmanned  Air  Vehicle  by  radar 
is  dependent  on  many  variables  including  range,  altitude,  and 
relative  orientation.  Given  a  radar  location  and  appropriate 
.model  for  the  likelihood  of  detection,  a  path  plan  can  be 
created  for  an  Unmanned  Air  Vehicle  which  constrains  the 
probability  of  detection.  In  this  paper  such  an  approach  Is 
taken  using  a  linearized  detection  model.  The  detection  model 
and  the  Unmanned  Air  Vehicle’s  dynamics  are  represented  as 
a  linear  program  subject  to  mixed  integer  constraints.  This 
mixed  integer  linear  program  (MILP)  Is  then  solved  with 
commercial  software  which  has  been  traditionally  used  by  the 
Operations  Research  community.  This  approach  searches  for 
all  feasible  solutions  and  produces  the  b^t  path  plan  based 
on  the  user  specified  parameters. 

1.  INTRODUCTION 

Optimal  path  planning  for  Unmanned  Air  Vehicles 
(UAVs)  is  a  difficult  task  due  to  the  intrinsically  non-convex 
nature  of  path  optimization.  Recently,  new  approaches  have 
been  used  which  utilize  techniques  traditionally  used  in  the 
field  of  Operations  Research.  By  using  linearized  dynamic 
models  and  mixed  integer  constraints,  researchers  have  used 
Mixed  Integer  Linear  Programming  (MILP)  for  different 
path  planning  applications.  [1]  [2]  [3]  MDLP  problem  for¬ 
mulations  can  be  solved  using  commercial  software  that 
employs  a  branch  and  bound  algorithm.  This  algorithm,  by 
looking  at  all  path  possibilities,  returns  the  optimal  solution, 
if  one  exists. 

Previous  efforts  [1]  have  concentrated  on  path  planning 
for  multiple  vehicles  and  multiple  goals.  Constraints  have 
included  time  to  waypoints  and  total  effort,  or  fuel.  In  this 
paper  the  MILP  approach  for  path  planning  is  expanded  to 
include  constraints  derived  not  only  from  the  relationship, 
but  also  the  interactions,  between  the  UAV  and  adversaries. 

Radar  guided  surface  to  air  missiles  (SAMs)  are  a  major 
threat  to  UAVs.  While  the  loss  of  a  UAV  does  not  result  in  a 
human  casualty,  it  can  jeopardize  the  mission  and  results  in 
the  loss  of  an  important  battlefield  resource.  Knowledge  of 
a  SAM’s  location  and  a  good  model  of  its  attack  capabilities 
can  be  used  to  create  a  UAV  path  plan  with  acceptable  risk. 
A  good  model  for  a  SAM’s  probability  of  detecting  a  UAV 
lakes  many  factors  into  consideration.  The  most  basic  is 
distance  between  the  SAM  and  the  UAV.  Another  important 
factor  in  determining  whether  a  SAM  site  detects  a  UAV 
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is  the  signature  presented  by  the  UAV  to  the  SAM’s  radar. 
Once  detected,  it  is  possible  for  a  UAV  to  avoid  destruction 
by  taking  advantage  of  the  possibility  of  lock  loss.  Lock  loss 
occurs  when  radar  detects  a  UAV  but  is  unable  to  continue 
detection,  or  lock  on,  long  enough  for  a  SAM  to  be  fired 
at  the  UAV. 

Non-convexities,  path  dependencies,  and  sharp  gradients, 
are  introduced  to  the  path  generation  problem  when  taking 
into  account  UAV  flight  dynamics,  the  probability  of  radar 
detection  based  on  UAV  state,  and  the  occurrence  of  lock 
loss.  Such  problems  have  been  previously  approached  using 
a  gradient  descent  based  method  which  can  result  in  local 
minima.  [4]  Discrete  approximations  to  continuous  shortest 
paths  have  also  been  investigated  in  relation  to  SAM 
avoidance.  [5]  The  approach  outlined  in  this  paper  uses  a 
linearized  model  of  ihc  probability  of  detection,  lock  loss 
phenomenon,  and  UAV  dynamics,  to  create  a  large  MILP 
formulation.  Once  solved,  this  MILP  formulation  results  in 
the  global  solution.  Examples  of  possible  path  plans,  as  a 
function  of  acceptable  probabilities  of  detection,  are  also 
presented. 

II.  General  Problem  Formulation 

The  model  is  based  on  the  Open  Experimental  Platform 
(OEP)  developed  by  Boeing[6].  Throughout  this  report,  it  is 
assumed  that  the  aircraft  maintains  a  fixed  altitude  and  that 
the  radar  is  on  the  ground.  The  model  is  presented  in  three 
parts:  the  vehicle  dynamics,  the  probability  of  detection 
model,  and  the  lock  loss  model. 

A.  Vehicle  Dynamics 

The  inputs  to  the  aircraft  model  include  measurements  of 
the  aircraft  position  and  velocity  as  well  as  the  destination 
waypoint  position.  A  constant  speed  for  travel  between  way- 
points  is  sAso  input.  The  output  for  the  aircraft  model  is  the 
aircraft  position  and  attitude  in  inertial  coordinates  (north, 
east,  up).  The  model  is  highly  simplified  and  does  not 
accurately  portray  the  true  physics  of  the  aircraft.  However, 
the  model  was  developed  for  research  in  mission  planning 
systems,  and  it  is  assumed  that  these  planning  controls  will 
be  designed  and  tested  in  conjunction  with  more  complex 
models  of  the  vehicles  and  their  flight  management  systems. 
The  model  assumptions  include  constant  altitude  (2D)  flight, 
instantaneous  changes  in  speed,  heading  and  bank  angle, 
and  simplified  turning  assumptions.  During  turns  the  bank 
angle  is  ±45°,  while  during  stead  level  flight,  the  bank  angle 
is  zero.  The  following  equations  are  used  for  the  aircraft 


state  for  ti  <t  <  U+i  where  U  denotes  when  the  UAV  is 
at  waypoint  i 

x{t)  =  x{ti)  +  U{ti)  COB{lp{ti))(t  -  u) 
y{t)  =  y{ti)  +  U{ti)  -  U) 

h{t)  =  h{ti) 

■4){t)  = 

m  ^  (1) 

where  x,  y,  and  h  are  the  aircraft  positions  along  the  north, 
east,  and  up  axes,  respectively  and  U  (t)  is  the  speed.  The 
heading  and  bank  angles  are  ip  and  ^  respectively.  The 
integer  variables,  Spt  and  Snty  are  defined  for  positive  (right) 
and  negative  (left)  turns  as  follows. 

Jpt  =  1  +  Tturn  and  V’(i<+i)  -  HU)  >  e 

Snt  =  l^t<ti  +  Ttum  and  ^(*<+1)  -  <  s 

where  e  is  a  smali  positive  constant.  The  heading, 
is  the  angle  between  the  nose  and  north  and  is  positive 
clockwise  about  the  up  axis  and  is  defined  below. 


ip{ti)  =  tan“^ 


The  steady  level  flight  turn  rate  equation  is  used  to  compute 
the  turn  time,  Tturn>  as  follows 


I  -  Hu)  I 

9tan(^) 

-mr 


To  compute  the  magnitude  of  the  turn  angle,  the  dot  product 
between  the  velocity  vector,  v{ti)  =  Vac(ft)[x]  -f  Vy(fi)[y], 
and  the  difference  between  the  destination  waypoint  and  the 
current  waypoint 


AR  =  (x{ti+t)  -  a;(*<))(£]  +  {y(<<+i)  -  y(*<))(yl 


is  used  as  follows 


cos(l  I)  = 


/  t^»(<i)(x(tf+i)  -  xjtj))  +  -  yjtj)) 

1’ 

_ 1 _ 

-  x{ti)f  +  (y(t<+i)  -  y{U)y 

For  the  sign  of  the  turn  angle,  the  cross  product  is  used. 
Specifically  if  y(ti^i)vx{ti)  —  a:(ii4.i)vy(ti)  >  0,  then 
the  change  in  turn  angle  is  negative,  ^(^*^.1)  —  ipiU)  = 
— l^(^«+i)  —  This  is  a  ’’left”  turn,  and  is  achieved 

by  rolling  the  aircraft  with  roll  angle,  <p  =  — 7r/4  rad, 
with  the  left  wing  down.  Similarly  if  y{ti^i)vx{ti)  - 
x{ti^i)vy(ti)  <  0,  then  the  change  in  turn  angle  is  positive, 
-  i’iU)  =  iV'C^t+i)  -  This  is  a  ’’right” 

turn,  and  is  achieved  by  rolling  the  aircraft  with  roll  angle, 
<p  =  7r/4  rad,  with  the  right  wing  down. 


Note  that  there  could  be  cases  where  there  is  no  turn. 
This  happens  if  the  turn  angle  is  less  than  a  small  amount 
specified  in  the  OEP.  Also,  there  could  be  cases  where  the 
turn  time  is  greater  than  or  equal  to  time  between  waypoints 
f*-f  1  —  ft-  In  this  case,  the  OEP  model  assumes  turning  for 
the  entire  segment. 


wmypoIntM 


AR  =  Wt^„)-x(t,))[x]+(y(t^„)-y(t,))[y] 
,)[xl+v,(0[y] 


-  v(i,)=v.(t,: 
^waypoint  t 


Hg.  1.  The  turn  angle  is  computed  from  the  velocity  vector  and  the 
vector  between  the  destination  and  current  waypoints. 


B,  Probability  of  Detection 

The  second  component  of  the  model  is  the  detection 
model.  This  model  computes  the  probability  of  detection 
of  a  UAV  by  an  opponent  SAM  radar.  The  inputs  of 
the  detection  model  are  the  outputs  of  the  aircraft  model, 
specifically,  the  aircraft  position  and  attitude.  The  positions 
of  opponent  radars  are  also  input  if  the  aircraft  is  within  an 
engagement  range  of  the  radar.  For  each  radar,  the  signature 
and  probability  of  detection  of  the  aircraft  is  included  in  the 
output  of  the  detection  model.  Signature  is  an  intermediary 
variable  that  is  related  to  radar  cross  section. 

For  each  radar,  a  vector  from  the  aircraft  to  the  radar 
is  transformed  to  aircraft  body  axis.  This  vector  is  then 
transformed  to  spherical  coordinates  for  azimuth  (Az), 
elevation  (El)  and  range  (R). 

The  signature  is  computed  with  a  table  look-up  as  func¬ 
tion  of  azimuth  and  elevation.  The  tables  used  in  this  paper 
are  shown  in  Table  1  and  Table  2.  Note  that  the  signature 
model  has  shaip  changes  between  azimuths  magnitudes 
between  30  and  31  degrees  and  also  between  159  and 
160  degrees.  When  the  nose  of  the  aircraft  points  directly 
to  the  radar  or  directly  away  from  the  radar,  and  the 
elevation  magnitudes  are  relatively  low,  there  is  a  region 
of  low  signature.  The  probability  of  detection  is  computed 
as  a  function  of  signature  and  range  using  another  table 
look-up.  For  the  probability  of  detection  model,  note  that 
probability  of  detection  increases  with  increasing  signature 
and  decreasing  range.  To  determine  probability  of  detection 
as  a  function  of  signature  and  range,  interpolation  is  used. 
The  maximum  and  minimum  probability  of  detection  are.99 
and.Ol,  respectively. 

In  the  OEP,  a  random  number,  uniformly  distributed 
between  0  and  1  is  generated.  If  this  number  is  greater 
than  or  equal  to  the  probability  of  detection,  the  UAV  is 
detected.  If  the  generated  number  is  less  than  the  probability 
of  detection,  the  UAV  is  not  detected.  We  assume  that  the 
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TABLE  I 

Signature  is  a  function  of  azimuth  and  elevation 
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TABLE n 

Probability  of  Detection  Table:  This  table  relates 

SIGNATURE,  RANGE  IN  KM,  AND  PROBABILITY  OF  DETECTION. 


Signature 


Probability  of  Detection 


PHH 

■EHI 

huh 

mumm 

mrsm 

HIHI 

g 

1 

■ 

mrsm 

mrmm 

m\mm 

■tVfM 

mnsurm 

mrsm 

mjmm 

miwm 

H££H 

mrsm 

HUJH 

miswm 

mnsm 

problem  is  deterministic,  and  that  a  detection  occurs  if  the 
probability  of  detection  is  greater  than  a  threshold  value. 

C  Lock  Loss 

A  mixed  logical  dynamical  (MLD)  representation  of  the 
lock  loss  model  is  presented  here  and  shown  in  Figure  2. 
Define  four  states:  disengage,  engage,  launch  and  damage. 
The  integer,  binary  variables  Sey  Sll,  ^launchy  ^damage 
used  to  signal  events  that  transition  the  opponent  SAM 
radar  system  from  one  state  to  another.  Whtn  Se  —  1, 
the  opponent  SAM  engages  the  UAV.  The  variables 
^launchi  ^damage  are  used  in  conjunction  with  timers  de¬ 
scribed  below  for  lock  loss,  launch  and  time-to-target.  When 
a  timer  surpasses  a  threshold  value,  a  change  in  state 
is  triggered.  When  the  lock  loss  timer  exceeds  threshold 
TiLy  variable  Sn  is  set  to  1  and  the  state  transitions 
to  disengage.  When  the  engage  timer  exceeds  Teng,  the 
variable  Siaunch  is  set  to  one,  indicating  that  the  SAM  is 
launching  a  missile  at  the  UAV.  When  the  time-to-taiget 
timer  exceeds  Ttargeu  the  variable  Sdamage  is  set  to  one, 
indicating  that  the  missile  has  had  enough  time  to  reach 
the  target.  These  conditions  are  detailed  by  the  following 
logical  equations  described  below: 


Se 

Sll 

^Launch 

^damage 


1  ^  engagement  conditions  met 

1  xll  >  Tll 

1  Xeng  ^  'I'eng 

1  Xiaunch  ^  '^target 

(3) 


The  variables  xn,  x^ng^  and  xiaunch  are  the  timer 
variables,  and  are  updated  using  the  equations  below. 


Fig.  2.  Lock  Loss  four  state  model. 


A  clock  is  defined  with  the  equation  below. 

i  =  0,l,2,... 

m  =  0  (4) 

Where  Ts  is  the  sample  time. 

For  the  lock  loss  timer,  if  there  is  no  detection,  the  lock 
loss  timer  is  incremented.  If  there  is  a  detection,  the  lock 
loss  timer  is  reset  to  zero.  The  lock  loss  threshold,  Tll>  ut 
(3)  is  specified  in  the  OEP  model. 

XLL{t{j  +  1))  -  XLL{t{j)){l  -  Se)  +  (1  -  Se)  (5) 

The  engage  timer  increments  when  the  system  is  engaged 
and  is  reset  to  zero  if  there  is  a  lock  loss  event. 

Xeng{^{j  +  1))  “  {^eng{^{j))  “F  ^e)(l  ~  (6) 

Similarly,  the  launch  timer  increments  when  a  launch  has 
occurred  and  is  reset  to  zero  if  there  is  a  lock  loss  event. 

Xlaunch{^{j  4"  1))  —  {xiaunch(j'{j))  “P  ^/o«ncA)(l  ^Ll)  (7) 


III.  MILP  Problem  Formulation 

In  this  section,  a  MILP  representation  of  the  problem 
is  presented.  The  first  subsection  pertains  to  developing 
a  vehicle  dynamics  model  that  is  suitable  for  the  MILP 
approach.  The  second  subsection  addresses  the  detection 
model,  and  the  third  pertains  to  the  lock  loss  model. 

In  order  to  represent  the  model  described  in  the  previous 
section  as  a  Mixed  Integer  Linear  Programming  (MILP) 
problem,  several  simplifications  are  made.  The  MILP  for¬ 
mulation,  like  the  model,  assumes  the  UAV  maintains 
constant  altitude.  In  addition,  the  bank  angle  of  the  UAV  is 
always  zero.  This  simplification  is  valid  if  the  turn  time  is 
small  compared  to  the  size  of  the  discrete  time  steps  used 
in  MILP. 
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i4.  Dynamics  Constraints 

The  total  flight  time  is  represented  throughout  the  for¬ 
mulation  as  T  discrete  time  steps,  which  is  a  parameter  the 
user  specifies  in  the  input  data  file.  The  UAV  is  given  T 
time  steps  to  reach  the  desired  end  point. 

There  are  T  values  for  the  UAV’s  x,  y  position.  These  T 
values  are  indexed  with  variable  fc.  At  any  individual  time 
step,  k,  the  UAV’s  position  at  that  time  step,  rrfA:]  and  y[k], 
is  calculated  to  be  the  sum  of  the  starting  x  and  y  position, 
called  xo  and  j/o»  plus  each  velocity  component,  Ua;[A:]  and 
Vy[k]  multiplied  by  the  length  of  time  of  each  time  step, 
which  is  a  constant  T,,  for  all  values  of  k  from  zero  to  the 
current  time  step.  Written  as  equations  this  gives: 


k 


x[k] 

=  lo  +  •  r,) 

3=0 

k 

(8) 

(9) 

Velocity  components,  Vx  and  Vy,  must  be  constrained  so 
that  the  magnitude  of  the  velocity  vector  is  not  greater  than 
the  maximum  velocity  of  the  UAV,  Vmox* 

Calculating  the  total  velocity,  using  Vx  and  Vy  results 
in  a  nonlinear  function,  due  to  the  square  root.  For  a 
MILP  representation  a  linearized  approach  must  be  taken 
to  constrain  Vx  and  Vy.  One  such  approach,  outlined  in  [1], 
is  to  constrain  Vx  and  Vy  separately  to  be  less  than  some 

Vm- 

I  Vx  |<  Vm  and  |  Vy  |<  Vm  (10) 

Plotted  on  a  Vx  versus  Vy  graph,  as  in  Fig.  3,  these 
constraints  create  a  box  that  is  inscribed  by  the  circle 
corresponding  to  all  feasible  Vx  and  Vy  values.  This  box 
is  a  poor  approximation  to  the  circle. 


vx 


Fig.  3.  Polygon  circle  approximation  for  M  =  4  and  M  =  8. 

To  achieve  a  better  approximation  to  the  circle,  more 
constraints  can  be  added. 


+u»[fclcos^^)  <Vn,a*cos(^)  (11) 

This  constraint  has  to  hold  for  all  integer  values  of  m  in 
the  set  from  1  to  M,  where  M  is  the  number  or  sides  of 
the  polygon  approximation  of  the  circle.  Since  the  polygon 
approximation  is  inscribed  by  the  circle,  Vmax  naust  be 
multiplied  by  the  cosine  factor.  V^thout  this  factor,  the  circle 
would  be  inscribed  by  the  polygon  and  there  would  be  Vx 
and  Vy  combinations  allowed  by  the  approximation,  which 
combine  to  produce  a  v  greater  than  Vmax^  This  is  easiest 
to  see  in  the  case  where  M  equals  four  since  there  is  a  large 
difference  between  the  box  inscribed  by  the  circle,  and  the 
box  that  inscribes  the  circle. 

The  above  constraints  only  ensure  that  v^ax  ts  not 
violated.  To  find  the  actual  velocity  at  any  time  step  fc, 
u[fc],  the  Vx  and  Vy  components  must  be  used.  Normally, 
calculating  the  magnitude  of  a  vector  from  its  components 
relies  on  the  nonlinear  square  root. 


vffc]  =  («*[*:])*  +  K[&])2  (12) 

A  linear  approximation  uses  three  constraints  and  several 
binary  variables,  by.  For  any  time  step  k  and  for  all  m  in 
the  range  of  1  to  M,  where  L  is  an  arbitrarily  laige  number: 

+Ui,[fc]co8^^^^  <  v[A:]  (13) 

+  Vy[k\ •  1-01  (14) 

>  v[k]  —  L(1  —  5v[fc,m]) 
and 

M 

5;;6„[fc,i]  =  i  ,  (15) 

i=i 

The  true  value  of  v[k]y  corresponding  to  the  UAV  velocity, 
will  satisfy  these  condition  for  only  one  value  of  m.  First 
(13)  constrains  v[k]  to  be  a  value  that  is  bigger  than  any 
possible  combination  of  Ux[fc]  and  Vy[k]y  for  any  value  m. 
Equation  (14)  also  constrains  u[fc]  to  be  smaller  than  any 
combination  multiplied  by  1.01,  or  else  the  corresponding 
binary  variable,  by[k]y  must  be  zero.  Wth  no  other  con¬ 
straints,  many  values  of  v[k]  could  satisfy  (13)  and  (14) 
by  having  all  zero  values  of  bv[k].  To  prevent  this,  (15) 
constrains  such  that  it  must  be  non-zero  for  one  m 
value.  This  m  corresponds  to  one  of  the  comers  of  the  M 
sided  polygon  appoximation.  When  bv[k]  does  equal  1,  (13) 
and  (14)  can  only  be  satisfied  with  a  t;[^]  which  is  a  close 
approximation  to  the  actual  velocity  value. 
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The  heading  value  is  important  in  determining  the  radar 
signature  of  the  aircraft.  The  heading  angle  is  tracked  with 
a  variable  which  will  later  be  constrained  to  ensure  that 
a  small  enough  cross  section  is  presented  to  any  radar  to 
avoid  detection.  Given  the  model  for  UAV  dynamics,  it  can 
be  assumed  that  the  UAV  always  faces  the  direction  it  is 
traveling.  Therefore  the  heading  is  related  to  the  velocity 
vector.  In  finding  v[k]  there  are  m  binary  variables,  by[k]^ 
for  each  time  step  k.  These  variable  all  equal  zero  except 
for  the  one  that  corresponds  to  the  m  which  relates  how 
Vx[k]  and  Vy[k]  combine  to  form  This  value  of  m  is 
directly  relat^  to  the  heading  angle  of  the  UAV  using  the 
following 

•  j  •  —  =  heading  angle  (16) 

J=i 

To  force  the  UAV  to  go  from  its  starting  position  to  the  final, 
or  desired  position,  more  constraints  and  binary  variables 
are  necessary.  Again  assuming  that  L  is  an  arbitrary  large 
value  we  can  write  the  following: 

x[k]-xf  <  L-(l-6/(A:])  (17) 

x[k]-Xf  >  -L-  {l-bf[k])  (18) 

y[k]-ys  <  i-(i-6#])  (19) 

yW  -  y/  >  -L  {\-  6/[A:])  (20) 

and 

T 

=  i  >  «'/b1e{0,i}  (21) 

i=i 

The  first  four  constraints  force  6/ [fc]  to  be  0  unless  the 
UAV  is  at  the  final  location,  Xf  and  j//,  in  which  case  bf[k] 
can  be  1.  The  final  constraint  says  that  when  the  binary 
variable  b  is  summed  over  the  entire  time  of  the  flight,  it 
must  equal  1.  These  constraints  force  the  UAV  to  the  final 
position. 

None  of  the  constraints  presented  so  far  force  the  UAV 
to  reach  the  final  position  as  quickly  as  possible.  To  do 
this  the  final  step  is  to  create  a  definition  for  what  makes  a 
path  optimal.  This  is  done  with  MILP  by  creating  a  metric 
that  characterizes  each  path  and  then  either  minimizing  or 
maximizing  that  metric.  Note  that  in  the  constraints  to  force 
the  UAV  to  the  final  position,  bf[k]  only  equals  1  when  the 
UAV  is  at  the  final  position  and  that  k  is  a  time  step  index. 
The  optimal  path  is  one  for  which  the  k  is  the  smallest. 
This  can  be  represented  by  defining  the  cost  as: 


Thus  the  optimal  path  is  the  one  that  satisfies  all  other 
constraints  and  has  the  lowest  cost. 


B,  SAM  Constraints 

Given  a  model  for  a  SAM’s  ability  to  detect  a  UAVj 
several  more  constraints  can  be  placed  on  the  UAV’s  path 
to  reduce  the  probability  of  detection  to  a  predetermined 
acceptable  level. 

As  described  in  the  model,  two  important  parameters 
in  determining  a  SAM’s  impact  on  the  UAV’s  path  is  the 
distance  between  the  two,  and  the  UAV’s  signature  to  the 
SAM’s  radar.  The  calculations  in  the  model  for  distance  and 
azimuth  are  nonlinear.  The  following  MILP  formulation  is 
an  equivalent  linear  approximation. 

EHstance  between  a  UAV  and  SAM  is  an  important  factor 
in  determining  the  probability  of  detection.  The  distance 
can  be  found,  and  then  constrained,  in  a  similar  fashion  as 
maximum  velocity,  using  constraints  and  binary  variables. 
For  any  time  step  k,  where  L  is  an  arbitrarily  large  number 
the  distance,  disf[A:],  is  constrained  as  follows, 

x[k]  -h  y[A:]  ^  dist[k]  (23) 

^x{A:]  sin^^  j  +  yj*]  j  j  •  1.01 .  (24) 

>  dist[k]  ~  L(1  —  6rf{fc,m]) 
and 

M 

EM*.y]  =  i  ,  Mfc.i]e{o,i}  (25) 

a;[A;]  and  y[k]  are  combined  so  that  they  are  less  than 
or  equal  to  disffA;]  yet  are  greater  than  or  equal  dist[k] 
when  multiplied  by  1.01.  This  must  hold  for  all  values  of 
m,  where  m  ranges  from  1  to  M.  The  binary  variables  are 
used  to  pick  out  the  one  true  value  of  distance  that  makes 
the  two  constraints  true. 

Another  factor  in  determining  whether  a  SAM  detects  a 
UAV  is  if  the  UAV  is  nose  in  or  nose  out  relative  to  the 
SAM.  Nose  in  refers  to  when  the  UAV’s  exposure  to  radar 
is  minimal.  Nose  out  refers  to  when  the  exposure  is  greater. 
This  is  determined  by  comparing  the  UAV’s  heading  angle 
and  the  line  of  sight  (LOS)  angle  between  the  SAM  and 
the  UAV.  The  LOS  angle  can  be  found  from  the  x,  y,  and 
distance  in  the  same  way  that  heading  was  found  from  Vx, 
Vy,  and  velocity.  Using  (25),  (26),  and  (27)  to  find  the  non¬ 
zero  value  of  bd[k]  the  line  of  sight  angle  is  given  by 

i]  =  Line  of  Sight  angle  (26) 

i=i  ^ 

With  the  line  of  sight  angle  and  the  previously  calculated 
heading  angle,  the  azimuth  angle  at  any  time  step  k,  can  be 
found  by  comparing  the  two.  The  line  of  sight  angle  and 
the  heading  angle  are  each  one  of  M  discrete  values,  where 
M  is  the  size  of  the  polygon  approximation.  A  static  M  by 
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M  table  can  be  used  to  find  the  azimuth  angle  as  a  function 
of  both  the  line  of  sight  angle  and  the  heading  angle.  To 
illustrate,  consider  if  the  UAV  has  a  heading  of  0°,  heading 
north,  the  azimuth  table  entry  would  be  180®  if  the  line  of 
sight  value  is  also  zero;  the  UAV  is  north  of  the  SAM  and 
heading  away.  However  if  the  line  of  sight  value  is  180®, 
than  the  azimuth  angle  would  be  0°,  the  UAV  is  south  of  the 
SAM  and  heading  directly  towards  it.  To  access  the  correct 
table  entry  binary  variables  bios  and  bh  are  used  as  indices. 

The  next  factor  to  consider  is  elevation  angle  between  the 
UAV  and  the  SAM.  Since  a  constant  altitude  is  assumed, 
the  elevation  angle  is  simply  a  function  of  range. 

dist  •  tan(elevation)  =  altitude  (27) 

However  tangent  is  a  non-linear  function.  Since  the  eleva¬ 
tion  table  given  in  the  model  for  SAM  detection  is  discrete, 
a  discrete  linear  approximation  can  be  used  for  the  tangent 
function  in  this  case.  For  any  time  step,  k,  an  index,  d, 
that  ranges  from  1  to  the  size  of  the  tangent  approximation 
array,  and  the  binary  variables  fceh  can  be  used  to  create  the 
following  constraints: 

alt  <  dist[k]  •  eltabie\d'  4- 1]  -h  M  •  (1  —  bei[ky  d\)  (28) 

alt  >  •  eftaWc(d]  M  •  (1  —  bei[ky  d])  (29) 

and 

£fre»(fc,i]  =  i  ,  (30) 

d=l 

^be,{k,d\-d--^  =  el[k]  (31) 

0=1 

Using  elevation,  azimuth,  and  acceptable  probability  of 
detection  set  by  the  user,  the  tables  given  in  section  2.2 
can  be  used  to  find  the  minimum  safe  distance  between 
the  UAV  and  the  SAM  radar.  In  the  MILP  formulation  the 
two  tables  given  in  section  2.2  are  combined  into  one  three 
dimensional  table.  First  the  elevation  angle  is  used  to  find 
the  correct  azimuth  versus  probably  of  detection  sub-table. 
Then  the  azimuth  variable  is  used  with  the  user  supplied, 
discrete,  probability  of  detection  threshold  variable,  to  look 
up  the  minimum  distance.  The  final  MILP  constraint  is  that 
the  UAV  distance  from  the  SAM  is  greater  than  this  look  up 
distance.  This  constraint  is  consistent  with  the  assumption 
that  if  a  path  violates  the  probability  of  detection  threshold, 
a  detection  is  assured.  This  simplification  is  necessary  to 
eliminate  the  random  element  of  SAM  radar  detection.  Due 
to  this  assumption,  for  any  path  to  be  valid  it  must  be  true 
that  for  any  given  time  the  distance  between  the  UAV  and 
the  SAM  is  greater  than  the  look  up  distance  based  on  all 
the  factors. 

C.  Lock-Loss  Constraint 

The  MILP  formulation  of  lock  loss  assumes  that  the  first 
time  the  UAV  is  detected,  the  lock  loss  timer  is  activated. 


During  the  lock  loss  time  if  the  vehicle  continues  to  be 
detected,  a  SAM  is  fired  and  the  UAV  is  destroyed.  If  the 
UAV  is  not  detected  for  as  long  as  the  lock  loss  time,  than 
the  engage  timer  is  reset  and  it  takes  another  detection  to 
start  it  again.  This  means  that  a  UAV  can  take  advantage  of 
lock  loss  by  flying  a  path  that  may  allow  it  to  be  detected, 
but  not  for  longer  than  the  lock  loss  time  limit.  In  MILP  this 
can  be  specified  using  binary  variables.  Suppose  the  lock 
loss  time  is  4  time  steps.  The  distance  constraint  is  changed 
to  include  6a,  binary  variables  for  the  lock  loss  such  that  at 
any  time  step  k  and  for  L  being  an  arbitrary  lai^e  number 

dist[k]  H-  (Ir  •  bii[k])  >  table  value  for  distance  (32) 

If  the  distance  at  time  k,  is  less  than  the  table  value  for 
distance,  this  constraint  can  only  be  true  if  bii[k]  is  1, 
otherwise  it  can  be  0.  To  ensure  lock  loss  the  final  constraint 
is  added,  assuming  the  lock  loss  time  is  4: 

3 

^6«Ifc+i]  =  3  ,  6H[fc  +  i]e{0,l}  (33) 

i=o 

Which  states  that  for  any  time,  k,  and  the  next  3  consecutive 
time  steps,  bu  must  be  equal  to  zero  once.  This  can  only 
happen  if  for  any  4  consecutive  time  steps,  the  distance 
between  the  UAV  and  the  SAM  is  greater  than  the  table 
value  at  least  once.  This  ensures  that  the  UAV  is  not  detected 
for  any  longer  than  the  lock  loss  time. 

IV.  Implementation 

The  AMPL  modeling  language  is  well  suited  for  rep¬ 
resenting  the  optimizations  presented  above,  A\^th  AMPL 
two  files  are  used,  the  model  file  details  the  various  con¬ 
straints.  Another  file  holds  the  model  data.  Thus  data  can 
be  changed,  such  as  starting  point  or  UAV  velocity,  without 
changing  the  constraints.  The  optimization  problem  coded 
in  the  two  AMPL  files  is  solved  with  ILOG  CPLEX  [7]. 
The  output  values  are  saved  to  a  third  file,  which  can  be 
opened  in  Matlab  for  data  analysis.  The  following  examples 
were  solved  on  a  1.8  GHz  Pentium  n  computer  with  512 
MB  of  RAM,  running  the  Windows  XP  operating  system. 

V.  Examples 

A.  Single  UAV  vs.  Single  SAM 

Using  the  data  in  the  tables  presented  in  Tables  I  and  n 
a  single  UAV  was  started  at  position  Xo  =  350  and  j/o  = 
125.  The  UAV  was  given  a  final  destination  of  Xf  =  340 
and  j//  =  0.  A  single  SAM  was  placed  at  the  origin  and 
an  acceptable  probability  of  detection  of  0.5  was  used.  The 
MILP  formulation  of  this  path  planning  problem  contains 
15363  binary  variables,  72  linear  variables  and  10359  linear 
constraints.  As  shown  in  Fig.  4,  the  UAV  does  not  fly 
directly  to  the  final  position,  but  instead  must  stay  far 
enough  away  from  the  SAM  until  it  reaches  a  position 
where  it  can  turn  towards  the  SAM  site  and  approach 
the  final  position.  The  circle  in  the  figure  represents  the 
minimum  distance  the  UAV  must  maintain  if  it  has  a  high 
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signature,  corresponding  to  a  nose  out  orientation.  The 
circle  appear  elliptical  due  to  the  unequal  scaling  of  x  and 
y.  The  UAV  only  turns  toward  the  destination  when  doing 
so  will  provide  a  small  cross  section  to  the  SAM’s  radar, 
allowing  it  to  get  closer  without  violating  the  acceptable 
probability  of  detection. 


B.  Effects  of  Probability  of  Detection 

By  varying  the  acceptable  probability  of  detection,  a 
more  direct  path  plan  can  be  generated.  Two  additional 
path  plans  were  calculated  using  all  the  same  parameters 
as  the  example  given  above,  except  that  the  probability 
of  detection  threshold  values  of  0.6  and  0.7  were  used. 
All  three  paths  are  plotted  in  Fig.  5,  which  clearly  shows 
that  when  a  higher  risk  of  det^tion  is  acceptable,  the 
UAV  can  use  a  more  direct  and  risky  path.  The  additional 
circles  correspond  to  the  smaller  minimum  distance  for  high 
signature  approaches.  The  smaller  the  probability  threshold, 
the  closer  the  UAV  can  be  to  the  SAM  at  the  origin  while 
presenting  it  with  a  high  signature. 

C.  Calculation  Time 

Each  path  given  in  Fig.  5  took  approximately  3  minutes 
to  compute.  This  time  was  needed  by  the  CPLEX  solver 
to  find  all  feasible  solutions.  Due  to  the  large  and  complex 
nature  of  the  MILP  formulation  there  are  many  branches 
for  the  solver  to  traverse  in  search  of  the  optimal  solution. 
While  not  fast  enough  for  realtime  applications,  it  should  be 
noted  that  this  calculation  time  is  greatly  influenced  by  the 
computational  capabilities  of  the  computer  used  to  solve 
the  MEvP.  Realtime  application  may  be  possible  with  a 
computer  optimized  to  solve  CPLEX  problems. 

D.  Effects  of  Lock  Loss 

Shown  in  Fig.  6  is  the  same  path  plan  from  Fig.  4,  with  a 
0.5  probability  of  detection.  Also  shown  is  another  path  that 
takes  advantage  of  the  lock  loss,  where  the  lock  loss  time  is 


Hg.  5.  UAV  paths,  SAM  at  origin  and  Probability  of  Detection  =  0.5, 
0.6,  and  0.7 


set  to  four  time  steps.  This  means  as  long  as  the  UAV  is  not 
detected  for  four  consecutive  time  steps,  the  UAV  will  not 
be  destroyed.  Due  to  the  extra  complexity  that  the  lock  loss 
modeling  adds  to  the  MILP  formulation  the  computation 
time  increased  to  approximately  3.5  mintues. 


Fig.  6.  UAV  paths,  lock  loss  and  non-lock  loss  with  Probability  of 
Detection  =  0.5,  SAM  at  origin. 


VI.  Conclusion 

As  the  use  of  UAVs  increases,  so  does  the  importance 
of  efficient  path  plans.  In  complex  situations,  such  as  the 
threat  of  a  SAM,  finding  a  feasible  path  plan  can  be  difficult, 
knowing  it  is  also  optimal,  even  more  so.  The  MBLP 
approach  outlined  in  this  paper  finds  the  most  efficient  path. 
Due  to  the  computation  time,  the  approach  outlined  above 
is  best  suited  as  a  top  level  path  planning  algorithm.  It 
may  also  be  possible  to  use  this  algorithm  in  a  receding 
horizon  fashion,  where  new  path  plans  are  generated  for 
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shorter  intermediary  waypoints  and  recalculated  as  new 
sensor  information  is  available. 

The  MDLP  approach  given  in  this  paper  does  not  require 
the  path  planning  to  be  computed  by  the  vehicle  itself.  One 
possible  scenario  could  entail  using  a  computational  facility 
to  solve  the  MILP  and  calculate  a  path  plan  that  is  then  sent 
over  a  communications  link  to  a  UAV. 
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